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1.  Introduction 


Extensions  of  the  original  constant-temperature  Dissipative  Particle  Dynamics  (DPD)  method 
have  been  developed;  in  particular,  methods  that  impose  constant-energy  (DPD-E)  ( 1 ,  2)  and 
constant-pressure  (DPD-P)  (3,  4)  conditions.  An  extension  not  currently  available  within  the 
DPD  framework  is  an  approach  that  imposes  constant-enthalpy  conditions  (DPD-H).  While 
conserving  total  enthalpy,  such  an  approach  would  allow  temperature  variations  within  the 
simulation  cell  under  constant-pressure  conditions.  Analogous  to  applications  of  constant- 
enthalpy  molecular  dynamics  (MD),  DPD-H  can  be  applied  to  pressure-dependent 
nonequilibrium  conditions,  e.g.,  where  mixing  becomes  prohibitive  as  the  simulation  cell 
pressure  increases  for  diffusion-limited  phenomena  (5).  As  part  of  the  work  presented  here,  we 
develop  a  DPD-H  method  by  combining  the  equations  of  motion  (EOM)  for  a  barostat  with  the 
EOM  for  the  DPD-E  method.  Both  deterministic  (Hoover)  and  stochastic  (Langevin)  barostats 
are  implemented  within  the  DPD-H  formulation,  where  a  barostat  temperature  is  defined  to 
satisfy  the  fluctuation-dissipation  theorem  (FDT)  for  the  Langevin  barostat. 

When  applying  any  DPD  variant,  numerical  integration  of  the  EOM  is  a  key  consideration 
because  the  stochastic  component  requires  special  attention.  Recent  work  by  our  group  (6-8)  has 
provided  a  comprehensive  description  of  numerical  integration  schemes  based  upon  the 
Shardlow-splitting  algorithm  (SSA)  (9)  for  the  isothermal,  isothermal-isobaric,  and  isoenergetic 
DPD  methods.  The  SSA  decomposes  the  EOM  into  differential  equations  that  correspond  to  the 
deterministic  dynamics,  and  elementary  stochastic  differential  equations  (SDEs)  that  correspond 
to  the  stochastic  dynamics.  In  the  original  SSA  formulation,  both  types  of  differential  equations 
are  integrated  via  the  velocity- Verlet  algorithm  (9),  where  the  stochastic  dynamics  are 
additionally  solved  in  an  implicit  manner  that  conserves  linear  momentum. 

In  this  work,  we  provide  an  SSA  formulation  for  the  DPD-H  method,  which  is  verified  by 
considering  both  the  standard  DPD  fluid  (pure  and  binary  mixtures)  (JO)  and  a  coarse-grain  solid 
model  (11).  For  both  models,  we  further  verify  the  DPD-H  variant  by  instantaneously  heating  a 
slab  of  particles  in  the  simulation  cell.  We  monitor  the  evolution  of  the  temperature,  pressure, 
and  density  as  the  system  approaches  an  equilibrated  state  while  maintaining  constant-enthalpy 
conditions.  For  completeness,  derivations  of  the  Fokker- Planck  equation  (FPE)  and  the  FDT  for 
the  DPD-H  variant  are  included  in  appendix  A. 
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2.  Formulations  of  DPD  at  Fixed  Pressure  and  Enthalpy  Using  Shardlow- 
Like  Splitting  Numerical  Discretization 


2.1  General  Formulation  of  DPD 

DPD  particles  are  defined  by  a  mass  m, ,  position  ri ,  and  momentum  p  .  The  particles  interact 
with  each  other  via  a  pair-wise  force  F„  that  is  written  as  the  sum  of  a  conservative  force  F.f , 

A  u  y 

dissipative  force  F(>D ,  and  random  force  F*  : 

F«=  Fs+F“  +  Fu-  (D 

F(  -  is  given  as  the  negative  derivative  of  a  coarse-grain  potential,  ufjG ,  expressed  as 


where  r;;  =  r,  -  r  is  the  separation  vector  between  particle  i  and  particle  j  ,  and  r.  =  |r;.  | .  The 
remaining  two  forces,  F.”  and  F;/R  ,  can  be  interpreted  as  a  means  to  compensate  for  the  degrees 
of  freedom  neglected  by  coarse-graining.  The  conservative  force  is  not  specified  by  the  DPD 
formulation  and  can  be  chosen  to  include  any  forces  that  are  appropriate  for  a  given  application, 
including  multibody  interactions  (e.g.,  11-13).  FyD  and  F*  are  defined  as 


and 


<4> 

where  y,(  and  a tj  are  the  friction  coefficient  and  noise  amplitude  between  particle  i  and  particle 

P  P; 

j  ,  respectively,  v,  =— - ,  and  Wr  are  independent  Wiener  processes  such  that  Wij  =  W ji . 

m.  m . 

’  J 

The  weight  functions  coD{r)  and  coR{r)  vanish  for  r  >  r, ,  where  r  is  the  cutoff  radius.  Note  that 
F^  is  completely  independent  of  Fj;D  and  F*  ,  while  F^  and  F*  are  not  independent  but  rather 

coupled  through  a  fluctuation-dissipation  relation.  This  coupling  arises  from  the  requirement  that 
in  the  thermodynamic  limit,  the  system  samples  the  corresponding  probability  distribution.  The 
necessary  conditions  can  be  derived  using  a  FPE,  which  are  derived  in  appendix  A  for  the  DPD- 
H  variant. 
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2.2  Constant-Enthalpy  DPD 

We  propose  a  constant-enthalpy  DPD  variant  based  upon  a  combination  of  the  EOM  for  a 
barostat  (7)  and  the  EOM  for  DPD-E  (8).  For  uniform  dilation,  the  combined  EOM  read  as 


dr,.  =^dl  +  —  rdf 


777; 


VE 
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dp,=2Xd(-r9 


00 
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-■yiJ 
r 

Vv  J 
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(5) 


where  Fs  =  dV(P -  P0)+  —  V^'  +  <tpW  ,  while  y.. ,  cr.. ,  ,  and  a.  ,  along  with  the 

/V,  ;  ///, 

weight  functions,  are  given  by  the  fluctuation-dissipation  relations  from  the  DPD-E  method  (8): 

al  =  2E^b©,; 

®D(r)=[®J!(r)]2 


E,  =  ^bE, 
(UD,?(r)=  [fo>AV/(r)]~ 


(6) 


where  the  relevant  temperature  is  ©3  =  — 

\ 

similar  to  coD{r)  and  coH(r),  respectively. 


1  1 

~0+~0 


,  and  coDl'{r)  and  aoRq{r )  can  be  chosen 


j  ) 
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Analogous  to  the  DPD-P  variant,  either  a  Hoover  barostat  ( yp  =  <jp  =  0)  or  a  Langevin  barostat 

can  be  invoked.  As  elaborated  in  appendix  A,  the  FDT  derived  from  the  FPE  for  the  Langevin 
barostat  parameters  yP  and  aP  read  as 

aP  =2  yPWekJbar,  (7) 


where  Tbar  is  the  Langevin  barostat  temperature,  which  is  an  additional  Langevin  barostat 
parameter  taken  together  with  yp ,  ap  and  W,: . 


The  system  temperature,  which  is  the  same  as  defined  in  DPD-E, 


T  = 


1  yP.P, 

3NkB  i  mi 


,  does  not 


necessarily  coincide  with  the  Langevin  barostat  temperature  since  the  DPD-H  variant 
corresponds  to  an  adiabatic  system.  More  generally,  the  equilibrium  state  (and  its  fluctuations)  of 
an  adiabatic  piston  in  contact  with  a  volume  reservoir  is  a  thermodynamically  ill-defined 
problem.  Effectively,  the  equilibrium  of  the  system  is  achieved  through  irreversible  processes 
(i.e.,  where  no  connection  with  an  external  reversible-work  device  such  as  a  Hoover  barostat 
exists)  that  depend  upon  the  dissipative  processes  occurring  in  both  the  heat  reservoir  and  in  the 
system  (14).  Consequently,  the  Langevin  barostat  parameters  given  in  equation  (7),  the  system 
compressibility,  and  the  system  viscosity  determine  the  amplitude  of  the  volume  fluctuations, 
and  thus  the  amplitude  of  the  system’s  energy  fluctuations.  When  Tbar  is  set  equal  to  the  system 


temperature,  the  equilibrium  state  (as  well  as  its  fluctuations)  are  identical  with  those  obtained  by 
the  Hoover  barostat.  The  effect  of  the  value  of  Tbar  on  the  system  properties  is  presented  in 


table  A-l  and  figure  A-l  of  appendix  A. 


Interestingly,  since  the  Hoover  barostat  corresponds  to  the  coupling  of  the  system  to  a  volume 
reservoir  through  a  reversible  mechanical  device,  the  volume  reservoir  is  neither  associated  with 
a  barostat  temperature  nor  does  it  produce  any  thermal  perturbations  on  the  system.  Rather,  the 
fluctuations  of  the  system  volume  itself  are  due  solely  to  thermal  perturbations  caused  by  the 
motions  of  the  particles,  and  thus  the  system  volume  fluctuations  are  independent  of  W£ . 


In  the  following,  we  demonstrate  that  the  EOM  for  the  Hoover  barostat  conserves  the  total 
enthalpy,  H  —  E  +  PQV  .  Analogous  to  constant-enthalpy  MD  (15),  DPD-H  with  the  Hoover 

barostat  should  also  conserve  the  quantity 


H'=H  + 


2  W 


=  e  +  p.y  +  ■ 


2W 


-Z 


p,  p, 

2/77, 


X  1  X  1  CG  ,  X  1  mech  ,  X  1  cond  ,  n  t  / 

Z.2AC  +23  +2Af  +poy- 


i  j>i 


2  W 


(8) 


4 


i.e.,  the  time  derivative  of  H'  should  be  zero.  First,  considering  the  time  derivative  of  E  only  in 
equation  (8)  gives 


d£_yP/  dp,.  yFc  dr,  |yd ur*  ,yd»r 

dr  ^/n,  dr  *  dr  ^  dr  4"  dr 


(9a) 


Substituting  the  corresponding  equations  given  in  equation  5  into  equation  9a,  and  using 

d  ucond 

y, — - —  =  0  (as  a  result  of  d uc°"d’l~j  =  -d uc°"d’I~j ,  cf.  equation  1  lc),  and  the  expression  for  the 
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Substituting  equation  9b  into 


d^=  d^+p  dF  |  pe  d pc 
dr  dr  0  dr  Ws  dr 


(9c) 


and  using  the  corresponding  equations  from  equation  5,  directly  leads  to  - =  0  ,  demonstrating 

dr 

that  the  EOM  in  equation  5  in  the  case  of  the  Hoover  barostat  conserves  total  system  enthalpy. 

Numerical  discretization  of  the  EOM  straightforwardly  follows  the  splitting  strategies  employed 
for  DPD-P  (7)  and  DPD-E  (5).  The  conservative  terms  of  the  deterministic  differential  equations 
are  identical  to  the  expressions  for  DPD-P  given  as 
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while  the  fluctuation-dissipation  terms  of  the  elementary  SDEs  are  identical  to  the  expressions 
for  DPD-E  given  as 
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As  before,  the  stochastic  flow  map  <f>At  is  approximated  from  the  first-order  splitting  algorithm 
given  by 
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(12) 


where  °  denotes  a  composition  of  operators. 

For  the  Shardlow-splitting  algorithm  velocity- Verlet  (SSA-VV)  approach,  updates  for  each 
term  are  performed  by 
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In  equations  13a  and  13b,  the  superscript  i  -  j  is  omitted  for  notational  simplicity.  Note  that  the 
total  system  energy  is  exactly  conserved  via  equation  13b.  In  equation  13a,  gq  -  -gqfi  is  a 

Gaussian  random  number  with  zero  mean  and  unit  variance,  chosen  independently  for  each  pair 
of  interacting  particles.  Finally,  </>%t  can  be  approximated  by  the  velocity- Verlet  algorithm 

given  by 
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where  L(o)  and  L(t  +  At)  are  the  lengths  of  the  cubic  simulation  box  at  t  =  0  and  t  +  At , 
respectively.  Next,  the  conservative  forces  at  t  +  At ,  (f f(t  +  A/  )}]V| ,  are  evaluated  and 
subsequently  used  in  the  second  part  of  the  algorithm,  which  requires  an  iterative  approach.  The 
iteration  starts  with  an  estimation  of  pe  at  t  +  At  using  p;  ']{l  +  At)  =  p£{t-  At)+  2A tFE(t) , 
followed  by  solving  the  set  of  equations 
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is  less  than  a  prescribed  tolerance,  which 


is  typically  less  than  10  6  .  In  equation  14b,  gp  is  a  Gaussian  random  number  with  zero  mean  and 
unit  variance  and  ( k )  is  the  iteration  index. 

A  practical  implementation  of  the  SSA-VV  for  DPD-H  follows.  The  EOM  and  the  numerical 
discretization  for  nonuniform  dilation  are  also  presented  in  appendix  B. 
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2.  Deterministic  Integration  #1.1 
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3.  Deterministic  Integration  #1.2  for  i  =  1  ,...,N 
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4.  Conservative  Force  Calculation:  {fc 


3.  Computational  Details 


The  SSA-VV  for  the  DPD-H  variant  was  tested  using  both  the  standard  DPD  fluid  (10)  and  a 
coarse-grain  solid  model  (11),  where  complete  details  of  the  conservative  forces  for  these  models 
are  given  in  appendix  C.  Both  a  pure  component  case  and  an  equimolar  binary  mixture  were 
tested  for  the  DPD  fluid  model.  System  sizes  for  the  DPD  fluids  and  coarse-grain  solid  were, 
respectively,  N  =  10125  and  13500.  For  these  simulations,  the  following  reduced  units  were 
used:  rc  and  r0  are  the  unit  of  length  for  the  DPD  fluid  and  coarse-grain  solid,  respectively;  the 
mass  of  a  DPD  particle  is  the  unit  of  mass;  and  the  unit  of  energy  is  kBTjni,  where  Tinj  is  the 
initial  system  temperature.  Using  these  reduced  units,  we  set  the  maximum  repulsion  between 
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particles  i  and  j  as  atj  =  25  for  the  pure  DPD  fluid,  and  as  atj  =  25  and  28  for  the  like  and  unlike 
i  -  j  interactions,  respectively,  for  the  binary  DPD  fluid.  Further,  we  set  the  noise  amplitude 
cr  =  3 ,  and  the  barostat  characteristic  time  zP  =  2 .  Prescriptions  for  the  choice  of  yp  (3, 16) 
suggest  that  the  value  should  be  between  Hzp  and  10/ rp ;  therefore,  we  set  yp  =  10/ rp  =  5  . 
Next,  we  assume  that  the  internal  degrees  of  freedom  are  purely  harmonic  and  express  the 
coarse-grain  particle  equation  of  state  as  u,  =  Cv  Jdi ,  with  heat  capacity  CV  i /kB=Cv/kB  =  60 
and  48  for  the  pure  DPD  fluid  and  coarse-grain  solid,  respectively.  Note  that  these  values 
correspond  to  coarse-graining  approximately  20  atoms  and,  separately,  16  atoms  into  a  DPD 
particle,  respectively  (17).  For  the  binary  DPD  fluid,  we  set  m2  =  10mj  and  Cv  2  =  10 Cv , ,  where 
Cv  |  /  kH  -  60 .  Finally,  following  Ripoll  et  al.  (18),  the  mesoscopic  thermal  conductivity  Ktj  is 
chosen  as 

K«=^}e<+e^  (15) 

where  k0  is  the  parameter  controlling  the  thermal  conductivity  of  the  DPD  particles,  which  are 
chosen  as  tc0  =  2.80  •  10^  for  the  pure  DPD  fluid  and  tcQ  =  1.52  ■  10-4  for  the  coarse-grain  solid. 
For  the  binary  fluid,  we  set  kb  u  -  2.80  •  lO^4 ,  kB22  =  2.80  •  10~6 ,  and  k0  12  =  xk0  22  .  Note  that 
when  6i  =  6 j  =  1 ,  these  particular  values  of  k0  and  Cv  give  Ktj  =  1 . 


4.  Results 


The  results  section  is  organized  as  follows.  The  validity  of  the  SSA  integration  algorithm  for  the 
DPD-H  variant  is  verified  by  considering  equilibrium  and  nonequilibrium  scenarios,  where 
results  are  given  in  subsections  1  and  2,  respectively.  Finally,  we  briefly  review  the  energy  drift 
associated  with  finite  integration  methods  and  propose  a  simple  strategy  to  minimize  these  drifts 
in  DPD-H  simulations. 

4.1  Test  Case  No.  1:  Equivalence  of  DPD  Variants 

As  a  first  test  of  the  SSA-VV  formulation  for  the  DPD-H  variant,  we  verify  that  it  converges  to 
the  same  equilibrium  properties  when  at  the  same  thermodynamic  conditions  as  other  DPD 
variants.  Starting  from  an  equilibrated  configuration  from  a  constant-temperature  DPD 
simulation,  DPD-P,  DPD-E,  and  DPD-H  simulations  were  performed  at  imposed  values  of  P  , 

E ,  and  H  ,  respectively,  where  these  values  were  determined  from  the  constant-temperature 
DPD  simulation. 

4.1.1  DPD  Fluid 

The  benchmark  systems  for  both  the  pure  and  binary  DPD  fluid  cases  are  taken  from  a  constant- 
temperature  DPD  simulation  performed  at  p  -  3  and  T  =  I .  and  run  for  trm  =  3000  and 
At  =  0.03 .  Further,  DPD-P  simulations  were  performed  using  both  the  Langevin  and  Hoover 
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barostats  with  an  imposed  pressure  determined  from  the  constant-temperature  DPD  simulation 
(P0  =  23.65  and  P0  -  24.79  for  the  pure  and  equimolar  binary  DPD  fluids,  respectively)  for 
trun  =  3000  and  At  -  0.03  .  The  following  quantities  were  evaluated  and  are  listed  in  table  1 
(pure  fluid)  and  table  2  (equimolar  binary  fluid):  virial  pressure  (Pir) ,  configurational  energy  per 
particle  (u) ,  kinetic  temperature  (Tkn\ ,  and  self-diffusion  coefficients  D  using  the  Einstein 
relation  (15). 

Table  1.  The  configurational  energy  per  particle  ( U  )  ,  the  kinetic  temperature  (Tkj„) , 
the  internal  temperature  (Tint)  » the  virial  pressure  ( Pvjr )  ,  the  particle  density 
(p  'j  ,  and  the  self-diffusion  coefficient  D  ,  determined  from  test  case  no.  1 
simulations  of  the  pure  DPD  fluid.  Tbar  is  the  Langevin  barostat  temperature 
and  (.')  denotes  an  ensemble  average,  where  numbers  in  parentheses  are 
uncertainties  calculated  from  block  averages. 


Variant 

(u) 

K) 

(Ktr) 

(P) 

D 

DPD 

P  =  3 

4.56(1) 

1.005(8) 

— 

23.65(8) 

— 

0.295(13) 

DPD-P 
Langevin 
P0=  23.65 

4.55(1) 

1.005(8) 

— 

23.59(8) 

2.997(8) 

0.294(14) 

DPD-P 

Hoover 

P0=  23.65 

4.55(1) 

1.004(8) 

— 

23.61(8) 

2.997(8) 

0.296(19) 

DPD-E 

P  =  3 

4.54(1) 

0.985(8) 

0.985(3) 

23.61(11) 

— 

0.293(6) 

DPD-H 

Langevin 

Tbar=  1-0 

P0  =  23.65 

4.54(2) 

0.985(8) 

0.985(3) 

23.65(11) 

3.002(11) 

0.294(7) 

DPD-H 

Hoover 

P0=  23.65 

4.54(2) 

0.986(8) 

0.985(3) 

23.65(9) 

3.002(13) 

0.292(8) 
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Table  2.  The  configurational  energy  per  particle  III  j  ,  the  kinetic  temperature  {j'kin  )  ,  the  internal 

temperature  (Tml)  ,  the  virial  pressure  (Pvir )  ,  the  particle  density  (p) ,  and  the  self-diffusion 
coefficient  D  ,  determined  from  test  case  no.  1  simulations  of  the  equimolar  binary  DPD 
fluid.  Tbar  is  the  Langevin  barostat  temperature  and  Q  denotes  an  ensemble  average,  where 
numbers  in  parentheses  are  uncertainties  calculated  from  block  averages. 


Variant 

(u) 

<r«> 

(p+) 

(P) 

A 

A 

DPD 

P  =  3 

4.76(1) 

1.005(8) 

— 

24.79(13) 

— 

0.177(13) 

0.165(13) 

DPD-P  Langevin 
PQ  =24.79 

4.76(2) 

1.004(8) 

— 

24.76(4) 

2.998(8) 

0.176(15) 

0.163(17) 

DPD-P  Hoover 

P0  =24.79 

4.76(1) 

1.004(8) 

— 

24.75(4) 

2.998(9) 

0.179(12) 

0.165(9) 

DPD-E 

P  =  3 

4.75(1) 

0.993(8) 

0.994(3) 

24.77(21) 

— 

0.174(9) 

0.161(10) 

DPD-H 

Langevin 

Tbar  =  1-0 

P0  =24.79 

4.75(2) 

0.992(8) 

0.994(3) 

24.78(13) 

3.001(9) 

0.174(8) 

0.160(9) 

DPD-H 

Hoover 

P0  =24.79 

4.75(2) 

0.993(8) 

0.994(3) 

24.78(19) 

3.001(13) 

0.173(8) 

0.160(10) 

To  validate  the  constant-enthalpy  SSA-VV  formulation,  DPD-H  simulations  were  performed  at 
conditions  taken  from  the  constant-temperature  DPD  simulation,  i.e.,  V  —  L'  —  3375  ( L  is  the 
box  length)  and  PQ  =  23.65  or  P0  —  24.79  for  the  pure  and  equimolar  binary  fluids,  respectively. 
DPD-H  simulations  were  carried  out  using  both  the  Langevin  and  Hoover  barostats,  where 
Tbar  =  1  was  used  for  the  Langevin  barostat.  The  final  configuration  of  the  constant-temperature 
DPD  simulation  is  used  to  determine  the  imposed  value  of  H ,  thus  it  is  also  used  as  the  starting 
configuration  for  the  DPD-H  simulation.  The  values  of  ut  were  initialized  by  setting  tr  =  Cv  Tmi, 
where  Tjni  =7  =  1,  and  were  carried  out  for  trun  =  3000  and  At  =  0.01 .  Analogous  to  isoenthalpic 
MD  simulations,  the  use  of  a  smaller  At ,  with  respect  to  constant-temperature  DPD  or  DPD-P 
simulations,  is  required  for  proper  conservation  of  H' .  Using  At  -  0.01 ,  we  observed  a  relative 
drift  in  H'  no  higher  than  1  •  10  4  .  (For  the  DPD  fluid  simulations,  reported  relative  drifts  refer  to 
an  average  of  relative  drifts  over  time  periods  of  1000.) 

Comparing  the  DPD-E  (5)  and  DPD-H  results  with  the  constant-temperature  DPD  and  DPD-P 
results  in  tables  1  and  2,  excellent  overall  agreement  is  found.  For  the  DPD-H  simulation,  the 

/  1  N  1  \  1 

internal  temperature,  (Tint)  =  /  — ^ —  )  was  also  evaluated,  where  the  values  of  (Tkin 'j  and  (Tmtj 

\N  i  ,  ()i  / 

agree  within  statistical  uncertainties.  For  the  pure  DPD  fluid,  these  values  are  approximately 
1.5%  lower  than  Tjnj  - 1 ;  however,  this  discrepancy  is  due  to  the  fundamental  differences 
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between  the  constant-temperature  DPD  and  DPD-E/DPD-H  methods.  Effectively,  the  two 
systems  are  different  since  the  imposed  temperature  for  the  constant-temperature  DPD  system 
should  be  equivalent  to  (Tkin ) ,  while  in  the  DPD-E/DPD-H  system  the  total  energy  initially  given 

to  the  system  is  dynamically  partitioned  among  the  kinetic  and  internal  energies,  yielding  a 
variation  in  the  equilibrium  temperature  with  respect  to  Tini  (19).  This  difference  is  of  0{kB  /  Cv) 


as  compared  with  unity,  while  an  additional  contribution  of  the  same  order  arises  from  the  “extra 
degree  of  freedom”  due  to  the  fluctuations  in  ui ,  since  the  relevant  macroscopic  temperature  is 


related  to 


'J_vl 


Hence,  for  a  pure  DPD  fluid  up  to  first  order  in  kB  /  Cv , 


Tkm '"'j  =  (Tint )  =  Tjnj(l  -kB/Cr)  =  0.983  (19),  in  agreement  with  the  simulated  values  of 

0.985  ±  0.003 .  For  the  equimolar  binary  fluid,  the  simulated  values  are  also  in  agreement  with 
the  estimate  (Tkm )  =  (Tinl)  =  0.994  .  Note  that  the  estimated  value  is  closer  to  Tini  =  1  due  to  the 


larger  value  of  Cv  2  that  reduces  the  (kB  / Cv)  contribution.  Due  to  these  lower  values  of  (Tkin 
and  ,  the  values  of  (p)  for  DPD-H  in  tables  1  and  2  differ  slightly  from  (p)  for  DPD-P 


Reproducing  equilibrium  averages  is  necessary  but  not  sufficient  proof  that  the  integration 
scheme  is  behaving  properly.  Hence,  as  a  further  demonstration  of  the  quality  of  the  SSA-VV 
and  the  proper  choice  of  At ,  for  the  pure  DPD  fluid,  we  calculated  probability  distributions  for 
Pj ,  Uj ,  and  V  for  constant-temperature  DPD  and  DPD-P  with  At  =  0.03 ,  and  for  DPD-H  with 
At  -  0.01 .  We  compared  the  probability  distribution  for  pt  with  the  corresponding  Maxwell- 
Boltzmann  distribution  (15),  while  the  probability  distributions  for  uj  and  V  were  compared 
with  those  obtained  with  a  very  small  At  -  0.001 ,  which  is  more  than  an  order  of  magnitude 
smaller  than  typical  values  of  At  used  and  is  thus  approximated  as  the  “exact”  result.  As  an 
example  of  results,  figure  1  presents  the  probability  distributions  of  pi ,  ui ,  and  V  for  DPD-H 
using  the  Langevin  barostat,  where  agreement  is  excellent. 
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Figure  1.  The  probability  distributions  of  (a)  p; ,  (b)  Uj ,  and  (c)  V  for  DPD-H 
with  a  Langevin  barostat.  Symbols  represent  simulation  results  for 
At  —  0.01 ,  and  lines  correspond  to  (a)  the  Maxwell-Boltzmann 
distribution  for  p:  and  (b,  c)  “exact”  results,  i.e.,  simulation  results 
where  At  —  0.001 . 
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4.1.2  Coarse-Grain  Solid 


A  validation  study  analogous  to  the  DPD  fluids  study  is  performed  for  a  coarse-grain  solid 
model,  where  a  recently  developed  nickel  model  is  considered  that  reasonably  reproduces  several 
measured  properties,  including  the  melting  temperature  (11).  For  a  benchmark  system,  a 
constant-temperature  DPD  simulation  is  performed  at  p  =  8260  kg/m3  and  T  =  1300K  for 

trun  =  1  ns  and  At  =  5  fs,  where  results  are  listed  in  table  3.  At  this  state  point,  the  atomistic 
Sutton-Chen  (SC)  model  of  nickel  predicts  a  pressure  of  approximately  0  bar  (11),  while  (Pvir) 
for  the  coarse-grain  solid  model  is  larger  than  0  bar.  Next,  starting  from  an  equilibrated 
configuration  from  a  constant-temperature  DPD  simulation,  nonuniform  dilation  DPD-P 
simulations  were  performed  using  both  the  Langevin  and  Hoover  barostats  at  P0  =  0  bar  for 

tmn  -  1  ns  and  At  =  5  fs,  where  results  are  also  shown  in  table  3. 

Table  3.  The  molar  configurational  energy  (ll')  ,  the  kinetic  temperature  {^kin)  ,  the  internal 
temperature  (Tjnr)  ,  the  virial  pressure  (^,,r)  ,  and  the  mass  density  {p) 
determined  from  test  case  no.  1  simulations  of  the  coarse-grain  solid  model  of 
nickel.  Tbar  is  the  Langevin  barostat  temperature  and  (}j  denotes  an  ensemble 
average,  where  numbers  in  parentheses  are  uncertainties  calculated  from  block 
averages. 


Variant 

(“) 

(kj/mol) 

{Pin) 

(K) 

(T,,) 

(K) 

<^*) 

(bar) 

(P) 

(kg/m3) 

DPD 

p  =  8260  kg/m3 

-508.44(11) 

1300.1(91) 

— 

5.91(94) 

— 

DPD-P 

Langevin 

P0  =  0  bar 

-508.43(14) 

1299.9(92) 

— 

-0.13(95) 

8259.3(73) 

DPD-P 

Hoover 

PQ  =0bar 

-508.44(14) 

1299.7(91) 

— 

0.06(99) 

8260.1(68) 

DPD-E 

p  =  8260  kg/m3 

-508.71(11) 

1274.8(88) 

1274.5(5) 

-150.49(98) 

— 

DPD-H 

Langevin 

Par  =  1 300  K 

P0  =0bar 

-508.86(14) 

1276.1(89) 

1275.7(11) 

-0.27(97) 

8273.2(73) 

DPD-H 

Hoover 

PQ  =0bar 

-508.87(13) 

1275.8(89) 

1275.6(5) 

0.07(986) 

8273.9(80) 
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Analogous  to  the  DPD  fluid  study,  DPD-H  simulations  were  performed  at  conditions  determined 
from  the  constant-temperature  DPD  simulation,  specifically,  P0  =  0  bar.  The  value  of  H  imposed 
in  the  DPD-H  simulation  was  determined  from  the  final  configuration  of  the  constant- 
temperature  DPD  simulation,  which  is  used  as  the  starting  configuration.  The  value  of  ut  was 
initialized  by  setting  ut  =  CvTjni,  where  Tjnj  =T  =  1300  K.  Nonuniform  dilation  DPD-H 
simulations  were  carried  out  using  both  the  Langevin  and  Hoover  barostats,  with  Thar  =  1300  K 
for  the  Langevin  barostat.  All  simulations  were  run  for  trun  -  1  ns  and  At  —  5  fs,  where  a  relative 
drift  in  H'  under  2  •  10  4  was  observed.  (For  the  coarse-grain  solid  simulations,  reported  relative 
drifts  refer  to  an  average  of  relative  drifts  over  time  periods  of  1  ns.)  Comparing  the  DPD-E  (5) 
and  DPD-H  results  with  the  constant-temperature  DPD  and  DPD-P  results  in  table  3,  excellent 
overall  agreement  is  found.  For  both  DPD-E  and  DPD-H,  the  values  of  (Tkin)  and  (Jlnt)  are  equal 
within  statistical  uncertainties.  These  values  are  approximately  2%  lower  than  Tini  =  1300  K  but 
agree  within  statistical  uncertainties  when  the  extra  degree  of  freedom  due  to  the  fluctuations  in 
ui  is  considered,  i.e.,  (Tkin )  -  (Tmt)  =  Tini(l-kB/ Cv)  =  1273  K  (19).  Furthermore,  due  to  these 
lower  values  of  (Tkin )  and  (Tnl') ,  the  values  of  {Pvir)  for  DPD-E  and  (p)  for  DPD-H  in  table  3 
differ  accordingly  from  (Pvir)  and  (p)  for  constant-temperature  DPD  and  DPD-P,  respectively. 

4.2  Test  Case  No.  2:  Heating  Response  in  DPD-H  Simulations 

As  a  second  test  case  to  verify  the  SSA-VV  formulation  for  the  DPD-H  variant,  a  nonequilibrium 
scenario  was  considered  for  both  the  DPD  fluids  and  the  coarse-grain  solid  model.  Starting  from 
a  final  configuration  of  a  constant-temperature  DPD  simulation,  a  slab  of  DPD  particles  in  the 
middle  of  the  simulation  box  was  instantaneously  heated  and  the  system  response  was  studied  at 
constant-  (P,H) conditions,  i.e.,  by  DPD-H  simulations. 

4.2.1  DPD  Fluid 

Analogous  to  test  case  no.  1,  the  final  configuration  from  the  constant-temperature  DPD 
simulation  (at  T  =  1  and  p  -3)  was  used  as  the  starting  configuration.  For  this  configuration,  a 
slab  of  particles  of  width  0.5L  in  the  middle  of  the  simulation  box  was  heated  by  assigning 
velocities  from  a  Maxwell-Boltzmann  distribution  corresponding  to  Theat ,  and  by  setting 
ul  =  CViTheat.  The  remaining  (nonheated)  particles  were  assigned  =  CViTini,  where  Ttm  =  T  =  1. 

As  a  test  of  the  DPD-H  variant  with  the  Langevin  and  Hoover  barostats,  a  simulation  was 
performed  using  Theat  =  10  and  P0  =  23.65  (pure  fluid)  or  P0  =  24.79  (equimolar  binary  fluid)  for 
Kun  ~  5000  and  At  =  0.005  (the  value  of  P0  corresponds  to  the  pressure  determined  from  the 
constant-temperature  DPD  simulation  at  T  =  1  and  p  =  3 ).  Tbar  -  Tini  was  used  for  the  Langevin 
barostat.  The  time  evolutions  of  Tkin  ,  Tint ,  Pvir ,  and  p  for  the  pure  DPD  fluid  are  shown  in 
figure  2,  where  a  relative  drift  of  6  •  10“5  in  H'  was  observed.  As  evident  in  figures  2a  and  2b, 
rapid  equalization  of  Tkin  and  Tiat  was  observed,  followed  by  increasing  Tkjn  and  7’nt  with  t ,  until 
reaching  equilibrium  at  t  >  600  .  Comparing  the  performance  of  the  Langevin  and  Hoover 
barostats,  the  Hoover  barostat  exhibits  substantial  oscillations  in  Pvir  and  p  at  the  beginning  of 
the  simulation  due  to  the  initial  heat  impulse,  after  which  the  oscillations  dampen  with  time. 
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Jakobsen  observed  similar  pressure  oscillations  (ringing)  for  the  DPD  fluid  using  a  Hoover 
barostat  (J).  Figure  2  shows  an  expansion  of  the  system  ( p  decreases  with  t )  to  compensate  for 
the  inputted  heat.  Early-time  behavior  is  displayed  in  the  insets  of  figures  2a  and  2b,  where  a 
sharp  decrease  in  Tkin  ,  due  to  the  interfacial  relaxation  of  the  cold  and  hot  regions,  and  the 
oscillating  response  of  Pvir  and  p  for  the  Hoover  barostat  are  shown  in  detail.  The  equimolar 
binary  DPD  fluid  exhibits  analogous  behavior,  where  the  values  of  (Tkin ) ,  (Tint) ,  (Pvir) ,  and  (p) 
determined  at  equilibrium  conditions  are  summarized  in  table  4.  Here  again  (Tkm )  =  (Tint)  is 
found,  along  with  equivalent  values  of  (Pvir)  and  (p)  for  the  Langevin  and  Hoover  barostats. 


17 


inset 


Figure  2.  Time  evolution  of  the  kinetic  temperature  Tkjn ,  internal  temperature  Tm  ,  virial  pressure 
Pvir ,  and  particle  density  p  for  a  DPD-H  simulation  of  the  pure  DPD  fluid  at  Pn=  23.65 , 
where  a  slab  of  particles  in  the  simulation  box  was  instantaneously  heated  by  Theat  =  10  at 
t  —  0  :  (a)  DPD-H  simulations  using  the  Langevin  barostat  with  a  barostat  temperature 
Thar  =  1.0,  and  (b)  DPD-H  simulations  with  a  Hoover  barostat.  Insets  display  early  time 
behavior  of  Tkm  ,  Tmt ,  p,  and  Pir . 
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Figure  2.  Time  evolution  of  the  kinetic  temperature  Tkjn ,  internal  temperature  Tm  ,  virial  pressure 
Pvir ,  and  particle  density  p  for  a  DPD-H  simulation  of  the  pure  DPD  fluid  at  P0  =  23.65 , 
where  a  slab  of  particles  in  the  simulation  box  was  instantaneously  heated  by  Theat  =  10  at 
t  —  0  :  (a)  DPD-H  simulations  using  the  Langevin  barostat  with  a  barostat  temperature 
Thar  =  1.0 ,  and  (b)  DPD-H  simulations  with  a  Hoover  barostat.  Insets  display  early  time 
behavior  of  Tkm  ,  Tm  ,  p  .  and  Prir . 


19 


Table  4.  The  kinetic  temperature  .  the  internal  temperature 

(7J  ,  the  virial  pressure  {Pvir )  ,  and  the  particle  density 

(p) ,  determined  for  Test  Case  #2  simulations  of  the 
pure  and  equimolar  binary  DPD  fluids.  Tbar  is  the 
Langevin  barostat  temperature  and  ()j  denotes  an 
ensemble  average,  where  numbers  in  parentheses  are 
uncertainties  calculated  from  block  averages. 

Pure  DPD  Fluid 


Variant 

K) 

K) 

(p*) 

(P) 

DPD-E 

P  =  3 

Thea,  =  10 

5.405(43) 

5.403(20) 

37.36(14) 

— 

DPD-H 

Langevin 

Tbar=  1-0 

P0  =23.65 

ThtM=\0 

5.374(41) 

5.372(20) 

23.65(11) 

2.247(8) 

DPD-H 

Hoover 

P0  =23.65 

ThtM=\0 

5.400(44) 

5.401(20) 

23.65(23) 

2.243(13) 
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Table  4.  The  kinetic  temperature  (Tjcin')  .  the  internal  temperature 
(T  t  ^  ,  the  virial  pressure  (Pvir )  ,  and  the  particle  density 
(p) ,  determined  for  Test  Case  #2  simulations  of  the 
pure  and  equimolar  binary  DPD  fluids.  Tbar  is  the 
Langevin  barostat  temperature  and  denotes  an 
ensemble  average,  where  numbers  in  parentheses  are 
uncertainties  calculated  from  block  averages  (continued). 

Binary  DPD  Fluid 


Variant 

(rj 

(p,) 

(P) 

DPD-E 

P  =  3 

Thea,  =  10 

5.340(44) 

5.350(16) 

38.36(14) 

— 

DPD-H 

Langevin 

^=1-0 

P0  =24.79 

Thea,=  10 

5.498(42) 

5.507(23) 

24.79(11) 

2.259(8) 

DPD-H 

Hoover 

P0  =24.79 

Theat=  10 

5.501(43) 

5.509(24) 

24.79(25) 

2.258(13) 

4.2.2  Coarse-Grain  Solid 

A  validation  study  analogous  to  the  DPD  fluid  study  is  carried  out  for  the  coarse-grain  solid 
model  of  nickel.  The  final  configuration  from  the  constant-temperature  DPD  simulation 
(at  T  - 1300  K  and  p  —  8260  kg/m3)  was  used  as  the  starting  configuration.  From  this  starting 
configuration,  a  slab  of  particles  of  width  0.5L  in  the  middle  of  the  simulation  box  was  heated  by 
assigning  velocities  from  a  Maxwell-Boltzmann  distribution  corresponding  to  Theat,  and  by 
setting  ui  =  CvTheat .  The  remaining  (nonheated)  particles  were  assigned  w  =  CvTini ,  where 
Tini  —  T  - 1300  K.  As  tests  of  the  DPD-H  variants,  simulations  were  performed  at  P0  =  0  bar, 
using  Theat  =  3000  K  for  trun  =  1  ns  and  At  =  5  fs.  For  the  DPD-H  simulation,  the  non-uniform 
dilation  approach  with  both  Langevin  and  Hoover  barostats  was  used;  Tbar  =  Tjnj  was  chosen  for 
the  Langevin  barostat.  Relative  drifts  of  1.8  •  10  4  in  H'  were  observed.  At  low  and  moderate 
pressures,  the  coarse-grain  solid  model  melts  between  1800  and  1850  K  (11).  As  a  result,  at  the 
end  of  the  DPD-H  runs,  the  particle  configuration  corresponds  to  a  liquid  state.  Figure  3  displays 
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the  time  evolution  of  Tkin  ,  7’nl ,  and  p  for  the  DPD-H  simulation  with  the  Langevin  barostat. 
Complete  melting  is  evidenced  by  reaching  a  plateau  in  the  time  evolution  of  p  for  the  DPD-H 
simulation,  where  complete  melting  occurs  at  approximately  0.15  ns.  Melting  at  constant-  ( V,E ) 
conditions  occurs  at  a  slightly  higher  temperature  (5)  due  to  the  pressure  build-up  within  the 
simulation  cell,  which  is  relieved  under  constant  ( P,H )  conditions.  The  values  of  (Tt;n ) ,  (Tint) , 
(Pvir) ,  and  (p)  determined  at  equilibrium  conditions  are  summarized  in  table  5,  where  again 
{Tkin)  =  {Tmt)  is  found. 


Figure  3.  Time  evolution  of  the  kinetic  temperature  Tkin  ,  internal  temperature  Tm  ,  and 
mass  density  p  for  a  DPD-H  simulation  using  a  Langevin  barostat  with  a 
barostat  temperature  Thar  =  1300  K  for  the  coarse-grain  solid  at  PQ  =  0  bar,  where 
a  slab  of  particles  in  the  simulation  box  was  instantaneously  heated  by 
Theat  =  3000  K  at  t  =  0 . 
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Table  5.  The  kinetic  temperature  (Tkjn  ) ,  the  internal  temperature  (7;nt)  ,  the 
virial  pressure  (^,,>)  ,  and  the  mass  density  l^p)  ,  determined  for 
Test  Case  #2  simulations  of  the  coarse-grain  solid  model  of  nickel. 
Tbar  is  the  Langevin  barostat  temperature  and  (}j  denotes  an 
ensemble  average,  where  numbers  in  parentheses  are  uncertainties 
calculated  from  block  averages. 


Variant 

(r,„) 

(K) 

(r„> 

(K) 

(n„) 

(bar) 

(P) 

(kg/m3) 

DPD-E 

p  =  8260  kg/m3 

7^=  3000  K 

2060.1(146) 

2060.0(28) 

8863.2(55) 

— 

DPD-H 

Langevin 

Ear  =  1  300  K 

P0  —  0  bar 

Ewat  =  3000  K 

2032.7(144) 

2032.1(9) 

0.0  (83) 

7339.3(89) 

DPD-H 

Hoover 

P0  =  0  bar 

7^=  3000  K 

2034.1(143) 

2033.8(8) 

0.3(76) 

7339.0(101) 

4.3  Conservation  of  Total  System  Enthalpy 

For  values  of  At  comparable  to  those  used  in  constant-temperature  DPD  and  DPD-P 
simulations,  we  observed  a  small  long-term  drift  in  H'  for  the  DPD-H  simulations.  For  example, 
for  the  values  of  At  used  for  the  DPD-H  simulations  in  this  work  (At  -  0.01  for  the  DPD  fluids 
and  At  —  5  fs  for  the  coarse-grain  solid),  the  small  relative  drift  produced  in  H'  was  typically  of 
order  10  4 .  When  At  was  decreased  by  an  order  of  magnitude,  the  relative  drift  in  H'  dropped 
to  order  10  7 .  A  typical  example  of  the  dependence  of  the  relative  drift  in  H'  on  At  for  the  DPD 
fluid  is  shown  in  figure  4.  The  values  of  other  properties  (not  shown  here),  such  as  the  kinetic 
and  internal  temperatures,  configurational  energy,  density  and  virial  pressure,  change  with  At  by 
less  than  0.5%.  This  behavior  is  comparable  with  microcanonical  or  isoenthalpic  MD  simulations 
when  the  velocity- Verlet  algorithm  is  used  for  the  integration  of  the  EOM  (15).  Since  the 
integration  of  the  fluctuation-dissipation  contribution  exactly  conserves  the  energy  (up  to 
machine  precision),  the  drift  is  caused  by  the  velocity- Verlet  algorithm  during  the  integration  of 
the  deterministic  contribution  in  the  DPD-E  or  DPD-H  EOM.  Similar  to  microcanonical  or 
isoenthalpic  MD,  a  long-term  drift  in  E  and  H'  is  thus  inevitable  in  the  SSA  for  DPD-E  and 
DPD-H. 
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Figure  4.  The  relative  drift  in  H '  as  a  function  of  the  integration  time  step  At  for 
DPD-H  simulations  with  the  SSA-VV. 


To  enforce  enthalpy  beyond  this  small  drift,  one  can  apply  the  following  numerical  procedure. 
After  each  time  step,  the  difference  between  the  current  H'  and  the  inputted  H'  is  calculated. 
This  difference  is  then  divided  by  the  number  of  particles  and  equally  subtracted  from  each  ui . 
This  is  a  useful  strategy  provided  that  the  drift  in  H'  has  a  mechanical  origin,  which  implies  that 
the  energy  drift  scales  as  kvT  .  Thus,  the  extra  energy  per  particle  subtracted  in  this  procedure  is 
very  small  compared  to  the  magnitude  of  ut ,  which  scales  as  C  T  .  In  this  work,  the  variation  of 
the  system  temperature  due  to  this  drift  was  found  to  be  negligible  and  the  dynamics  unaffected. 
This  strategy  was  applied  to  all  test  cases  for  DPD-H  with  a  Hoover  barostat,  where  no  variation 
in  the  results  was  observed.  Note  that  H'  is  not  a  fixed  quantity  for  DPD-H  with  a  Langevin 
barostat  (rather  it  fluctuates  about  an  average  of  H');  therefore,  this  error  suppression  scheme 
can  only  be  applied  to  DPD-H  with  a  Hoover  barostat. 
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5.  Conclusion 


An  isoenthalpic,  isobaric  DPD  method  was  developed  by  combining  the  equations  of  motion 
from  both  the  DPD-P  and  DPD-E  methods.  Both  Hoover  and  Langevin  barostats  were 
implemented,  where  an  additional  barostat  parameter,  the  barostat  temperature,  was  defined  for 
the  DPD-H  variant  with  the  Langevin  barostat.  A  comprehensive  description  of  a  numerical 
integration  scheme  based  upon  the  Shardlow- splitting  algorithm  was  presented  for  the  DPD-H 
approach,  which  was  found  to  be  a  readily  extendable  and  accurate  integration  scheme. 

The  equivalence  of  the  DPD-H  variant  was  verified  using  both  a  standard  DPD  fluid  model  and  a 
coarse-grain  solid  model,  where  thermodynamic  quantities  as  well  as  probability  distributions 
were  considered.  The  integration  algorithm  was  further  verified  by  considering  equilibrium  and 
nonequilibrium  simulation  scenarios.  Finally,  a  discussion  of  the  inevitable  small,  long-term  drift 
in  H'  associated  with  finite  integration  methods  was  given,  where  we  propose  a  simple  strategy 
to  minimize  the  effect  of  this  drift  in  DPD-H  simulations. 
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Appendix  A.  Fokker-Planck  Equation  (FPE)  and  Fluctuation-Dissipation 

Theorem  (FDT) 


The  FPE  corresponding  to  the  equation  of  motion  (EOM)  given  by  equation  5  of  the  technical 
report  is 


Lcp  +  LDp  +  Lcondp  +  LLBp , 
ot 


(A-l) 


where  p  =  p(r,p,  V,  p£,u;t).  The  conservative  operator  Lc  is  given  by 
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The  operator  LD  representing  the  effects  of  the  dissipative  and  random  forces  is  given  by 
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The  operator  Lcond  associated  with  the  effects  of  the  mesoscopic  heat  transfer  between  particles 
is  given  by 
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a l  =  0  to  ensure  that  the  EOM  contain  no  spurious  drift. 


The  operator  Llb  representing  the  Langevin  barostat  terms  in  the  EOM  is  given  by 
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E  J 


Analogous  to  the  constant-energy  (DPD-E)  variant,  an  implicit  heat  reservoir  is  not  present  in 
constant-enthalpy  Dissipative  Particle  Dynamics  (DPD-H),  nor  is  the  total  system  enthalpy 
specified  in  the  FPE  (A-l).  Therefore,  to  simplify  the  derivation  of  the  FDT  for  DPD-H,  we 
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choose  the  isothermal-isobaric  ensemble.  Under  the  steady-state  condition,  peq  =  peq(r,p,V,p£,u ) 
corresponds  to  the  probability  density. 1,2,3 


peq(r,p,V,p£,u)  =  ^exp 


X— +  +  P0V  +  ^ 

^  2/77 ;  11  2  W 
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*Br 


(A-6) 


Analogous  to  the  previous  cases  6-8,  Lcpeq  =  0  for  the  equilibrium  distribution.  The  FDT  then 
follows  from  the  requirements  that  Lnp"‘  -  0 ,  Lcondpeq  =  0  and  LLBpeq  -  0 ,  which  lead  to 

equations  (6)  and  (7)  given  in  the  report.  As  mentioned  in  the  report,  while  the  FDT  relations 
relevant  to  the  dynamics  of  the  particles  depend  only  on  particle  variables,  the  FDT  relevant  to 
the  Langevin  barostat  depends  on  the  heat  reservoir  temperature.  For  constant-pressure  DPD,  this 
temperature  is  simply  the  system  temperature;  however,  in  constant-enthalpy  DPD,  the  heat 
reservoir  temperature  (i.e.,  barostat  temperature)  becomes  an  additional  Langevin  barostat 
parameter  along  with  yp ,  ap,  and  W£  ,  which  can  be  written  as 

cr2P=2yPW£kJbar,  (A-7) 


where  Tbar  is  the  Langevin  barostat  temperature.  Interestingly,  a  Hoover  barostat  is  a  purely 
mechanical  barostat  with  no  dissipation  (nor  an  associated  barostat  temperature);  hence,  the 
volume  fluctuations  are  only  a  result  of  the  inherent  thermal  perturbations  from  the  motions  of 
the  particles.  In  contrast,  a  Langevin  barostat  incorporates  dissipation  via  the  friction  and  random 
terms  in  Fe  (see  equation  5  of  the  report);  thus  the  volume  fluctuations  of  the  system  will 
additionally  depend  on  the  value  of  Tbar .  (As  an  illustrative  analogy,  one  can  consider  the 
Hoover  barostat  as  a  piston  in  a  vacuum,  while  the  Langevin  barostat  can  be  viewed  as  a  piston 
in  a  viscous  fluid,  whose  temperature-dependent  viscosity  affects  the  piston  fluctuations.) 

To  demonstrate  the  effects  of  Thar  on  the  system  properties,  we  carried  out  DPD-H  simulations 
using  the  Langevin  barostat  with  Tbar  =  {0.5,1.0,1.5,2.0}for  the  DPD  fluid.  We  also  performed 
DPD-H  simulations  using  the  Hoover  barostat,  where  all  simulation  results  are  summarized  in 
table  A-l.  The  table  clearly  shows  that  the  values  of  the  volume  fluctuation  SV2  and  the 
associated  isentropic  compressibility  depend  on  Tbar ,  while  thermodynamic  properties  such  as 
the  configurational  energy  per  particle,  kinetic  temperature,  internal  temperature,  virial  pressure, 
and  particle  density  are,  as  expected,  not  affected  by  Tbar .  Additionally,  figure  A-l  displays  the 


Jakobsen,  A.  F.  Constant-Pressure  and  Constant-Surface  Tension  Simulations  in  Dissipative  Particle  Dynamics.  /.  Chem. 
Phys.  2005,  122,  124901. 
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^Mackie,  A.  D.;  Bonet  Avalos,  J.;  Navas,  V.  Dissipative  Particle  Dynamics  With  Energy  Conservation:  Modelling  of  Heat 
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Quigley,  D.;  Probert,  M.  I.  J.  Langevin  Dynamics  in  Constant  Pressure  Extended  Systems.  J.  Chem.  Phys.  2004, 120,  1 1432. 
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volume  probability  distributions  for  the  Langevin  barostat  at  different  Tbar  and  for  the  Hoover 
barostat,  which  further  elaborates  the  effects  of  Tbar  on  SV2 .  In  table  A-l  and  figure  A-l,  the 
volume  fluctuations  corresponding  to  the  Langevin  barostat  with  Tbar  =  1.0  and  those 
corresponding  to  the  Hoover  barostat  are  nearly  equivalent  since  Tbar  =  1.0  is  close  to  the 
equilibrium  system  temperature  of  approximately  0.986.  If  Tbar  was  set  to  the  value  of  the 
equilibrium  system  temperature,  the  volume  fluctuations  given  by  the  Langevin  barostat  and 
those  given  by  the  Hoover  barostat  would  match  precisely. 

Table  A-l.  The  configurational  energy  per  particle  (u^  ,  the  kinetic  temperature  (^fen  )  ,  the  internal 
temperature  (7;nt)  ,  the  virial  pressure  \Pvir)  .  the  particle  density  (pj  ,  the  isentropic 
compressibility  Ks  ,  and  the  volume  fluctuation  SV  2  for  the  DPD  fluid  at  various  values 
of  the  Langevin  barostat  temperature,  Tbar.  denotes  an  ensemble  average,  where 
numbers  in  parentheses  are  uncertainties  calculated  from  block  averages.  These 
simulations  were  all  started  from  the  same  final  configuration  of  a  constant-temperature 
DPD  simulation  at  T  =  1  and  particle  density  p  =  3 .  The  imposed  pressure  P0  =  23.65 
corresponds  to  (Pvir )  from  the  constant-temperature  DPD  simulation.  Since  Cv  /  kB  =  60 
was  used  in  the  coarse-grain  particle  equation-of-state,  equilibrium  temperatures  should  be 
1  —  kB/Cv=  0.983 ;  Cv  is  the  heat  capacity. 
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Figure  A-l.  The  probability  distributions  of  (a)  Pt ,  (b)  LL  ,  and  (c)  V  for 

DPD-H  with  a  Langevin  barostat.  Symbols  represent  simulation 
results  for  At  —  0.01 ,  and  lines  correspond  to  (a)  the  Maxwell- 
Boltzmann  distribution  for  p{  and  (b,  c)  “exact”  results,  i.e., 
simulation  results  where  At  —  0.001 . 


33 


Intentionally  lelt  blank. 


34 


Appendix  B.  Constant-Enthalpy  Conditions  (DPD-H)  for  Nonuniform 

Dilation 


Following  the  derivation  of  the  uniform  dilation  constant-enthalpy  conditions  (DPD-H)  variant, 
the  combined  equations  of  motion  for  non-uniform  dilation  using  a  Langevin  barostat  and 
constant-energy  conditions  (DPD-E)  are  given  as 
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(u/)2dt  (i=l,...,iV) 


where  d  is  the  dimensionality  of  the  system  and  the  variables  ,  pea,  W,. ,  V  ,  Fea,  and  Pafj 
are  the  same  as  in  the  non-uniform  dilation  DPD-P  variant,  except  oy>  =  2 yPWEkBTbar,  and  Tbar  is 
the  Langevin  barostat  temperature.  Variables  pertaining  to  u"ie,h  and  u""ul arc  the  same  as  those 
in  the  uniform  dilation  DPD-H  variant. 


B.l  Numerical  Discretization 

Applying  a  numerical  integration  splitting  strategy  similar  to  the  uniform  dilation  DPD-H 
variant,  the  deterministic  differential  equations  and  the  elementary  constant-energy  (SDEs) 
corresponding  to  equation  (B-l)  are  the  following.  The  conservative  terms  are  the  same  as  those 
for  the  non-uniform  dilation  constant  pressure  conditions  (DPD-P)  variant,1  while  the 
fluctuation-dissipation  terms  are  identical  to  the  expressions  given  in  equation  1 1  of  the  report. 


*Bonet  Avalos,  J.;  Mackie,  A.  D.  Dissipative  Particle  Dynamics  With  Energy  Conservation.  Europhys.  Lett.  1997,  40  (2),  141. 
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As  stated  in  the  report,  the  stochastic  flow  map  is  approximated  by  equation  12  of  the  report. 

For  each  fluctuation-dissipation  term  ,  momenta  and  internal  energies  are  updated  by  using 
the  same  expressions  as  in  the  uniform  dilation  DPD-H  variant  for  the  constant  energy 

(SSA-VV)  scheme  (equations  13a  and  13b  of  the  report).  can  be  treated  in  identical  fashion 
to  the  nonuniform  dilation  DPD-P  variant. 
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For  the  models  considered  in  this  work,  the  details  of  the  conservative  forces  expressed  in 
equation  2  of  the  main  text  are  the  following.  °  for  the  pure  and  binary  constant-temperature 

Dissipative  Particle  Dynamics  (DPD)  fluids  is  given  by 

(C-l) 

where  is  the  maximum  repulsion  between  particle  i  and  particle  j  . 


For  the  coarse-grain  solid  model,  which  has  a  face-centered-cubic  (f.c.c.)  lattice  structure, 
particles  interact  through  a  shifted-force  Sutton-Chen  embedded  potential  (SC)  given  as: 


=  £ 


'kiyr-'i a 

Z  j*i  J 


(C-2) 


where 


(C-3) 


s  and  r0  are  the  energy  and  length  parameters,  respectively,  n  and  m  are  positive  integers 
(n  >  m  to  satisfy  elastic  stability  of  the  crystal),  and  c  is  a  dimensionless  parameter.  Although 
effectively  this  is  a  many-body  potential,  the  force  on  each  particle  can  be  written  as  a  sum  of 
pair-wise  contributions.  The  coarse-grain  solid  model  used  here  approximates  nickel  (Ni),  where 
one  constant-temperature  Dissipative  Particle  Dynamics  (DPD)  particle  was  chosen  to  represent 
four  f.c.c.  unit  cells,  i.e.,  16  Ni  atoms.  SC  potential  parameters  were  determined  by  fitting  to 
various  0-K  properties  and  the  melting  temperature  at  zero  pressure,1  where  the  following  values 
were  found:  s/kB  =  225  K,  r0  =  8.8698  A,  c  =  39.4314 ,  m  =  6 ,  and  n  =  9 .  Further  details  for 
determining  SC  parameters  based  upon  such  a  procedure  can  be  found  elsewhere  and  references 
therein. 1 


'Brennan,  J.  K.;  Lfsal,  M.  Proceedings  of  the  14th  International  Detonation  Symposium,  Coeur  d’Alene,  ID,  11-16  April 
2010;  Office  of  Naval  Research,  2010,  p  1451. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


DPD 

DPD-E 

DPD-H 

EOM 

f.c.c. 

FDT 

FPE 

MD 

SC 

SSA 


constant-temperature  Dissipative  Particle  Dynamics 

constant-energy  Dissipative  Particle  Dynamics 

constant-enthalpy  Dissipative  Particle  Dynamics 

equations  of  motion 

face-centered-cubic 

fluctuation-dissipation  theorem 

Fokker-Planck  equation 

molecular  dynamics 

Sutton-Chen  embedded  potential 

Shardlow- splitting  algorithm 
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